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O ABSTRACT 

Numerical calculations of merging black hole binaries indicate that asymmetric emission of gravita- 
tional radiation can kick the merged black hole at up to thousands of km/s, and a number of systems 
have been observed recently whose properties are consistent with an active galactic nucleus containing 
a supermassive black hole moving with substantial velocity with respect to its broader accretion disk. 
We study here the effect of an impulsive kick delivered to a black hole on the dynamical evolution 
of its accretion disk using a smoothed particle hydrodynamics code, focusing attention on the role 
played by the kick angle with respect to the orbital angular momentum vector of the pre-kicked 
r \ disk. We find that for more vertical kicks, for which the angle between the kick and the normal 

vector to the disk 6 < 30°, a gap remains present in the inner disk, in accordance with the prediction 
1-Q from an analytic collisionless Keplerian disk model, while for more oblique kicks with > 45°, 

Mh matter rapidly accretes toward the black hole. There is a systematic trend for higher potential 

Q luminosities for more oblique kick angles for a given black hole mass, disk mass and kick veloc- 

^ ity, and we find large amplitude oscillations in time in the case of a kick oriented 60° from the vertical. 

C^ Subject headings: Black hole physics. Accretion, accretion disks. Hydrodynamics 
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1. INTRODUCTION 



In the past few years, a combination of surprising numerical studies and astronomical observations have indicated 

^H that asymmetric momentum losses in gravitational radiation from the mergers of binary black holes (BH) may produce 

C^^ "kicks" of up to several thousand km/s. For the mergers of supermassive black holes (SMBH) at the centers of galaxies, 

^^ the kicks may be large enough to eject the remnant BH out of the galactic center if not out of the galaxy entirely. 

f^^ Kicked BHs may have already been observed indirectly as active galactic nuclei (AGN) with different components at 

^^ different redshifts; broad-line regions that are thought to originate near the BH itself will remain bound to the kicked 

^H BH and acquire its new line-of-sight velocity, while narrowline systems that are produced further away will become 

^H unbound and remain behind. 

K^ While kicks from the mergers of unequal-mass, non-spinning BHs have long been predicted by post-Newtonian 

• ^-H calculations (Fitchett & Detweiler 1984), the ability to evolve black holes in fully general relativistic simulations 

S^ has considerably expanded our view. Indeed, it is the merger of equal-mass, spinning BH that produce the largest 

^ kicks calculated to date, with maximum values of up to 4000km/s possible for configurations with carefully chosen 

^ alignments (Campanelli et al. 2007a, b; Gonzalez et al. 2007). While speeds this large represent only a small fraction of 

the likely merging BH parameter space, Monte Carlo simulations of merging BHs with arbitrary spins of dimensionless 

magnitude SjM^ = 0.97 find that the mean kick for BH systems with mass ratios uniformly distributed between 

and 1 is 630km/s, with more than 20% of the kicks larger than lOOOkm/s (Lousto et al. 2010). Such kicks have the 

potential to unbind the remnant from smaller galaxies, or displace the BH and any bound gas for larger galaxies. In 

roughly half of all major mergers between comparable-mass SMBHs, the merged SMBH should remain displaced for 

30 Myr outside the central torus of material that would power an AGN (Komossa & Merritt 2008) . 

Somewhat more recently, several AGN have been observed to contain broadline emission systems that appear to 
have very different redshifts than the narrowline emission systems. In the first of these, SDSSJ092712. 65+294344.0, a 
blueshift of 2650km/s is observed for the broadline systems relative to narrow lines (Komossa et al. 2008). Although 
several different physical models have been proposed to account for the observations, including broadline emission from 
the smaller SMBH within a binary surrounded by a disk (Dotti et al. 2009a; Bogdanovic et al. 2009), a pair of SMBHs at 
the respective centers of interacting galaxies (Heckman et al. 2009), and even spatial coincidence (Shields et al. 2009a), 
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the kicked disk model remains entirely plausible. Since then, several more systems with similar velocity offsets between 
different AGN emission regions have been discovered, including SDSS J105041. 35+345631. 3 [3500km/s; (Shields et al. 
2009b)], CXOC J100043.1+020637 [1200km/s; (Civano et al. 2010)], and E1821+643 [2100 km/s; (Robinson et al. 
2010)]. 

Given a number of recent results that suggest that binary SMBH spins should align with the angular momentum 
of the accretion disk prior to merger, these large kick velocities are somewhat unexpected. Indeed, aligned spins 
can produce kicks of several hundred km/s (Campanelli et al. 2007b) but not the "superkicks" with recoil velocities 
> lOOOkm/s. Bogdanovic et al. (2007) suggest that the accretion torques in gas-rich ("wet") mergers should suffice 
to align the SMBH spins with the accretion disk prior to merger. In Dotti et al. (2009b), high resolution simulations 
of SMBH binaries whose orbits counter-rotate with regard to a surrounding disk indicate that they should undergo 
an angular momentum flip long before merger. With typical spin-orbit misalignments of no more than 10° — 30° 
depending on the parameters of the disk, in particular its temperature, they find that the recoil kick during mergers 
should have a median value of 70km/s, with superkicks an exceedingly rare event (Dotti et al. 2010). 

While some of the observed candidate recoil velocities are so large that they fall well out into the tail of the kick 
velocity probability distribution function, it is difficult to constrain exactly how many systems indicate potential recoils 
given the challenges in clearly distinguishing multiple velocity components within a single AGN. Given this difficulty, 
we attempt here to study potential electromagnetic (EM) signatures that would originate from post-merger disks. 

The qualitative details of pre-merger evolution have been studied by other groups, and a relatively coherent picture 
emerges. In general, each SMBH may be surrounded by a circum-BH accretion disk extending out to a distance 
substantially smaller than the binary separation, since orbits near the outer edge are unstable to perturbations from 
the other SMBH. A circumbinary disk may be present as well, extending inward to a distance of roughly a few times 
that of the binary separation, with a gap appearing between the circumbinary disks and the inner disks. Previous 
hydrodynamical studies have shown that the inner edge of the circumbinary disk is driven into a high-eccentricity 
configuration that precesses slowly, while a 2-armed spiral density wave is formed extending out to larger radii (Mac- 
fadyen & Milosavljevic 2008). Meanwhile, the inner disks will be fed by mass transfer from the circumbinary disk 
(Hayasaki et al. 2007; Hayasaki et al. 2008), as the m = 2 azimuthal gravitational perturbation induces an elongation 
in the outer disk. For circular orbits, the mass transfer rate is relatively constant, while for elliptic binaries the mass 
transfer takes on a more periodic character. Finally, once the binary begins its gravitational radiation-driven plunge, 
the binary decouples form the outer disk, and mass transfer basically ceases. 

Among the first predictions of the observational consequences of a post-kick accretion disk are those in Loeb (2007), 
who find that off-center quasars could be observed for up to 10^ years after a SMBH merger. Schnittman & Krolik 
(2008) find that the inner edge of the circumbinary disk is likely to occur at roughly lOOOM, in typical relativistic 
units where G = c = 1, with the value only weakly dependent on the typical disk parameters like the assumed a- 
parameter. Based on semi-analytic models of the post-merger accretion disk, they find the potential for observable 
infrared emission lasting hundreds of thousands of years, leading to the prediction that several such sources might be 
present today, though they would be difficult to disentangle from other AGN sources. On much shorter timescales 
immediately prior to the merger, the dissipation of gravitational radiation energy through spacetime metric-induced 
shearing of the disk can also lead to enhanced emission in the optically thin components of the disk (Kocsis & Loeb 
2008). 

One of the first studies of post-kick disk dynamics was performed in Shields & Bonning (2008), who analyzed the 
approximate physical scales characterizing the post-merger disk and concluded that the total energy available to be 
dissipated by shocks is roughly ^M^v^-^j^, where M5 is the mass of the portion of the disk that remained bound and 
^kick is the SMBH kick velocity. They estimated from a simple analytical model that an excess luminosity of ~ 10^^ 
erg/s would be observable for roughly 3000 years, with a characteristic observed temperature of ~ lO^K assuming 
'^kick = lOOOkm/s. This model is checked by means of a collisionless N-body simulation of a disk around a kicked 
SMBH, which found rough agreement with the analytical model and confirmed the prediction of a visible soft X-ray 
flare that would last for a few thousand years after the initial kick. 

Using a 2-d version of the FLASH code, Corrales et al. (2010) are able to study the response of a thin disk if the 
kick is directed in the equatorial plane. They find that the characteristic response is a one-armed spiral shock wave, 
capable of producing total luminosities up to ~ 10% of the Eddington luminosity on a timescale of months to years. 
The relativistic decrease of the total SMBH mass, attributable to the energy carried off in gravitational waves and 
roughly 5% of the original total, leads to a decrease in the luminosity of approximately 15% but does not provide 
a clear signature that can be disentangled from the other global processes occurring in the disk. Nevertheless, with 
future X-ray instruments such as the Square Kilometer Array, impulsive changes to a disk may be observable in the jet 
emission (O'Neill et al. 2009). Even in the case where the BH kick is very small, the secular, as opposed to dynamical, 
filling of the inner region of the disk should produce an afterglow that could exceed the Eddington luminosity if the 
accretion rate is sufficiently large prior to the binary decoupling and merger (Shapiro 2010). 

The most direct comparison to the calculations we present here is found in Rossi et al. (2010). Using an analytic 
treatment of disks with power-law density profiles, they construct a model for disk evolutions in which, immediately 
after the kick, fiuid elements are assumed to circularize at radii determined by their specific angular momentum, with 
the resulting energy gain by shocks released as EM radiation. Their work establishes that the primary energy reservoir 
for kicked disks is potential energy that can be released by elliptical orbits, not the relative kinetic energy of the 
kick itself nor the impulse sent through the disk by instantaneous mass loss to gravitational radiation by the central 
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SMBH during the merger. They also perform detailed 2-d Eulerian and 3-d SPH simulations of post-kick accretion 
disks, although there are some important differences between the latter simulations and most of our SPH studies, as 
we discuss throughout the paper below. Among these, they assume that disks are isothermal, with all shock heating 
immediately radiated away, whereas we evolve the energy equation and allow the fluid to heat. They also include an 
"accretion radius", Racc^ such that particles that approach closer than that distance to the black hole are accreted 
and removed from the simulation domain (G. Lodato, private communication), whereas all particles remain in our 
simulation throughout the duration of a run. As such, their results and ours bracket the range of possible heating 
scenarios. 

The 2-d calculations by Rossi et al. (2010) of razor-thin disks using the ZEUS code indicate that vertical kicks 
(that is, in the direction of the disk angular momentum vector) lead to modulated emission, unlike all other kick 
angles they consider. This result is confirmed by their 3-d SPH simulations. In-plane kicks develop a clear spiral- wave 
structure with accretion streams forming as the simulation progressed, but a smooth luminosity profile. Their 3-d 
SPH simulations with oblique (less vertical) kicks indicate essentially a 2-phase model for the observed EM emission. 
Immediately after the kick, the majority of the luminosity may be attributed to the innermost region of the disk 
dissipating the kinetic energy it acquired during the kick, while at later times, after roughly one standard timescale of 
the disk (t = 1 in the notation of Eq. 1 below) and thereafter, the dominant dissipation mode is potential energy from 
infalling material on elliptical orbits. Post-kick disks are found to be rather compact, extending out to roughly the 
"bound radius" (f = 1 according to the notation of our Eq. 2) with steep density dropoffs at larger radii. Luminosities 
are generally highest for in-plane kicks, with roughly a factor of four difference in peak luminosity between largest 
peak luminosity (in-plane kick) and the smallest (vertical kick), with peaks occurring later for more oblique kicks. 
Relativistic BH mass decrease was found to be unimportant at large radii, and potentially important only in the 
vicinity of the BH (out to a few hundred Schwarzschild radii) where the effect of the kick is merely perturbative 
compared to the nearly relativistic Keplerian velocities. 

Numerical relativity groups have also considered the hydrodynamics of matter around both binary SMBHs and 
kicked BHs, typically at much smaller size scales. In Bode et al. (2010), flows of gas around a binary SMBH system 
are considered at scales roughly 10^ smaller than typical Newtonian calculations, spanning scales roughly lAU across, 
rather than ~parsec scales. They find that EM emission is dominated by variability created by Doppler beaming of 
the SMBHs as they shock the gas, leading to an EM signal that demonstrates the same periodicity as the gravitational 
wave signal, with corresponding peaks in the timing of the maximum emission for each. In a follow-up work. Bode 
and collaborators (Bode et al. 2011) predict that observable EM emission from near the SMBHs is much more likely 
to arise in a hot accretion flow, in which a flare would be seen coincident with the merger. In Megevand et al. (2009), 
the effect of a kicked BH moving through the equatorial plane of an accretion torus is considered using a fully general 
relativistic Eulerian hydrodynamics code. In their simulations, the newly-merged SMBH is surrounded by a torus 
extending out to 50M (50AU for a IO^Mq SMBH), and the overall timescale studied is approximately 10,000M (in 
relativistic units with G = c = 1), or about two months. They find that a kick in the direction of the equatorial 
plane of the torus produces the strongest shock in the system and therefore the strongest EM emission, consistent with 
studies that examine disks on substantially larger scales. By ray-tracing their simulations in a post-processing step 
(Anderson et al. 2010), they confirm that simple Bremsstrahlung luminosity estimates yield a qualitatively accurate 
picture of the disk luminosity for high-energy radiation, while their torus is optically thick to low-energy emission. 

Numerical calculations of vacuum EM fields surrounding an SMBH merger indicate that they could contribute to 
periodic emission (Palenzuela et al. 2009, 2010b) but are likely to be too small in amplitude and at the wrong frequencies 
(~ 10~^Hz) to be observed directly (Mosta et al. 2010). Such mergers could produce observable levels of Poynting flux 
in jets (Palenzuela et al. 2010a), however, through a binary analogue of the Blandford-Znajek mechanism (Blandford 
& Znajek 1977), which seems especially effective for spinning BHs (Neilsen et al. 2010). Calculations of non-Keplerian 
accretion disks in general relativity indicate that the spiral wave structures seen in Newtonian simulations could exist 
in relativistic models with small disks when the BH kick is sufficiently small, but that larger kicks disrupt the spiral 
pattern, as could dissipative processes such as magnetic stress or radiative cooling (Zanotti et al. 2010). The inferred 
emission due to synchrotron emission from a relativistic disk is considered by the Illinois group (Farris et al. 2010, 
2011), who find that emission could peak at a luminosity of ~ 10^^ erg/s, a few orders of magnitude brighter than the 
corresponding bremsstrahlung luminosity and potentially observable by either WFIRST or LSST. 

The outline of the paper is as follows. In Section 2, we introduce the physical scales that define the kicked disk 
problem, and discuss the basic dynamics of disks around a kicked black hole when hydrodynamic effects are ignored. In 
Sec. 3, we describe the SPH code used to perform 3-dimensional simulations, including versions that seek to replicate 
previous results and those that use new techniques to bracket the range of physical physical predictions. Results from 
these simulations, focusing on the hydrodynamic and thermodynamic evolution of the disk, are reported in Sec. 4. 
Finally, in Sec. 5, we lay out consequences of our results and plans to extend them in the future. 

2. PHYSICAL SCALES AND QUALITATIVE EXPECTATIONS 

2.1. Physical scales 

Throughout this paper, we use scaled units, denoted by hats, under the assumption that G = Mbh = '^kick = 1. A 
single unit of time, for example, is thus GMbuv^^^^. Choosing reference values Mbh = 10^ Mq and Vkick = 10^ cm/s. 



Ponce, Faber, & Lombardi 



the resulting time and distance scales are defined as 
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Because the disk is evolved without self-gravity, the disk mass scale is formally independent of the BH mass that sets 
the rest of the physical scales. 

For a disk of total mass mdisk, we may define a set of quantities marked with tildes by choosing a reference value 
^disk = lO^M©. The physical scales for energy, its time derivative, volume density, and surface density are then given 
respectively by 
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Assuming the fluid is a fully ionized ideal gas with mass fractions X = 0.7, Y = 0.28, the mean molecular mass is 
/i = 2/(1 + 3X + 0.5F) = 0.617, and the characteristic temperature scales like 

ks V 10 cm/s y 

where rup is the mass of a proton, and ks is the Boltzmann constant. The optical depth for Thomson scattering, 
r = hie^^ where Kq ~ 0.2(1 + X) = 0.34 cm^/g is appropriate for non-relativistic plasma, is given by 
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2.2. Qualitative expectations 

The simplest model for a disk around a merging SMBH binary consists of a 2-dimensional, inflnitely thin disk with a 
perfectly Keplerian rotational profile. If we assume that Newtonian gravity applies, we have a nearly scale-free system, 
where only the mass of the SMBH binary contains units. In what follows, we work in dimensionless units such that 
G = Mbh = 1, where Mbh is the total mass of the SMBH binary, assumed for the moment to be equal to the total 
mass of the merged SMBH that will be formed in the merger and immediately kicked with velocity I'kick at an angle 
relative to the angular momentum of the 2-d disk. We also choose '^kick = 1 to set the overall scaling for the remaining 
quantities we consider, denoting all quantities in these units by hats. The Keplerian rotational velocity, for instance, 
is given by vq = t~^I^ . Note that our unit system results in the speed of light not being set to unity. Note that our 
conventions are very similar to those found in Rossi et al. (2010), except for the definition of the kick angle: we define 
Q to be the angle away from vertical, whereas they define d to be an the angle between the kick and the initial disk 
plane. 

In what follows, we assume that the disk orbits in the x — y plane, with angular momentum in the positive z-direction. 
We define to be the azimuthal angle in this plane, measured counterclockwise from the x axis. The BH kick falls in 
the X — z plane, and we assume that the SMBH proceeds to move off with constant velocity, feeling no accelerations 
from either the disk or the galactic potential on the timescales of interest. 

The relative velocity and specific energy of points in the disk with respect to the newly kicked BH is given by 
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Fig. 1. — Particle periastron values after the BH kick as a function of radius for azimuthal angles ^ = —60°, — 90° and selected kick 
angles Q. Assuming the disk has an inner edge at r = fo prior to the kick, the minimuni value of fp decreases as B increases from 0°, a 
vertical kick, to ^ = 90°, an in-plane kick, with particularly rapid changes for large values of Q. 

Accordingly, the critical radius f^ as a function of azimuthal angle ^ inside of which the matter is bound to the kicked 
BH and unbound outside (cf. Fig. 1 of Rossi et al. (2010)) is given by 



n 



1 



(4) 



For both bound and unbound particles, we may calculate the full trajectory of the particles assuming they remain 
collisionless by solving the two-body problem with respect to the kicked black hole. Of particular importance for our 
later discussion in the post-kick periastron radius, fp, shown in Fig. 1, for longitudes ^ = —60° and —90° where the 
values can be quite small, indicating significant potential mass flow toward the BH. Values are systematically smaller 
for kick angles closest to the original disk plane {d = 90°), and, if we impose an inner edge on the original disk, there is 
a minimum periastron for all post-kick particles that decreases dramatically as the kick angle becomes more in-plane. 
This leads to a testable prediction for 3-d dynamical simulations: a "gap" should be present for more vertical kick 
angles, disappearing only once collisions in the inner region of the post-kick disk facilitate the angular momentum loss 
required to feed flow toward the SMBH, as noted for in-plane and vertical kicks in Rossi et al. (2010). For vertical 
kicks, we expect that it should take several times the orbital period at the inner edge of a disk before any signiflcant 
amount of matter is present near the SMBH. For in-plane kicks, the fllling of the inner disk should be much more 
rapid. 

3. SMOOTHED PARTICLE HYDRODYNAMICS FOR KICKED ACCRETION DISKS 

3.1. Initial data 

Following the methods outlined in Rubio-Herrera & Lee (2005a, b) and similar works, we construct semi-analytic 
models of accretion disks in hydrodynamic equilibrium to use as initial conditions before laying down particles using 
a Monte Carlo technique. To do so, we flrst assume a that the orbital velocity is independent of the height within the 



6 Ponce, Faber, & Lombardi 

disk and varies only with cylindrical radius. Integrating the force equation for a system in stationary equilibrium, 



P 

and assuming the pressure P = kp^ ^ where 7 = 5/3 is the adiabatic index of the fluid, M the mass of the central BH, 
r and Vc are the spherical and cylindrical radii, Tc is the unit vector associated with the cylindrical radius, and livc) 
is the chosen radial specific angular momentum profile, we find 

T] = -— = + / — :^drc - K, (6) 

7 — 1 p r J rg 



where K is an integration constant. For the Keplerian profile livc) = \fGMT~c^ the integral on the right hand side of 
eq. (6) evaluates to —GM/vc^ and the disk can exist only in the z = plane (where r = re). 

To get a disk with finite extent in the radial and vertical directions, we choose a non-Keplerian rotation profile. 
Here, we assume a power-law form, noting that it must satisfy l{r) oc r'^ with tz < 0.5 to yield a compactly bounded 
configuration. We define rotation parameters through the relation 

%^rfr, = -< (7) 
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and find that the top edge of the disk, where /? = 0, yields the condition that 



Assuming (hatted) units such that G = M = 1, the inner and outer edges of the disk for a sub-Keplerian rotation 
profile {n < 0.5 and thus a < — 1) are given by the two real roots of the equation z{r) = 0, or 

cf"+^+i^f = l. (9) 

To fix the inner and outer radii at f^ and tq respectively, we determine c{a) and K{a) as follows. Defining Ri = r^~^^ 
and Ro = f^^"*^, we find 
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For the case a = —2, corresponding to a constant l{rc)^ the solution is easy to state explicitly: noting that Rq = l/r^ 
and Ri = 1/ri, we find 
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In general, the easiest method to achieve a specific disk height Zmax is to vary a and check numerically where the 
disk reaches its maximum height, iterating until the correct value is achieved. Relatively thin disks for which the radial 
extent is significantly greater than the vertical tend to be nearly Keplerian, with a = — 1 — eQ,, where the eQ, is positive 
and < 1. 

For all the runs shown below, we chose initial parameters f^ = 0.1, fo = 2.0, and Zmax = 0.2 for our initial disk, 
resulting in an SPH discretization and rotation curve we show in Fig. 2. This corresponds to a choice of parameters 
c = 0.9584, a = —1.017, and K = 2.651 x 10~^. Once we specify the adiabatic index 7 = 5/3, we are left with a free 
parameter in the adiabatic constant k = P/ p^ . Varying k has the effect of rescaling the density (see eq. 6), and thus 
allows us to adjust the overall disk mass while leaving a uniform specific entropy. 

3.2. Hydrodynamical Evolution 

The code used to perform our calculations combines the parallelization found in the Star Crash SPH code with a 
number of refinements described in Lombardi et al. (2006) and Gaburov et al. (2010), which we summarize in an 
appendix. Among these, we implement variable smoothing lengths for all of our particles following the formalism 
described in Springel & Hernquist (2002); Monaghan & Price (2001); Price & Monaghan (2007), which can be derived 
consistently from a particle-based Lagrangian. Gravitational softening for particles near the BH is implemented in a 
self-consistent manner, similar to the method used in Price & Monaghan (2007), which is critical to avoid individual 
and collective particle numerical instabilities for matter in the vicinity of the BH, as we discuss below. 
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Fig. 2. — Top: The rotational profile of the disk, expressed as a ratio of the actual tangential velocity to the the Keplerian velocity, 
shown as a function of radius. Bottom: Initial disk configuration, projected onto the re — z plane. It has an inner radius f^ = 0.1, outer 
radius fo = 2.0, and maximum height Zmax = 0.2. Finite disks require super-Keplerian rotation on the inner regions and sub-Keplerian 
rotation outside so that the centripetal acceleration respects the additional pressure force component. The density maximum occurs at 
f = 0.23, where vt = Vkep- 

Artificial viscosity is included in all runs using a standard "Balsara switch" (Balsara 1995). In most of our runs, we 
implement the standard thermodynamic energy equation and allow the fluid to heat. Such cases provide a minimum 
estimate for the luminosity that can be produced by a kicked disk, since thermal pressure will damp continued shocking 
within the disk by lowering the ambient densities. Since this ignores the effects of radiative cooling losses, we also 
consider the more common approximation made in numerical simulations (see, e.g. Corrales et al. (2010); Rossi 
et al. (2010), that all shock-generated energy is immediately radiated away, resulting in an adiabatic evolution (or 
isothermal, depnding on the initial hydrodynamic model). Details of the implementation are discussed in the Appendix 
(see Eq. A3). While it is relatively straightforward to calculate the luminosity in such a model, the accretion of large 
numbers of particles very near the BH is an inevitable consequence, and further code modifications are required to 
allow numerical evolutions, particularly the introduction of an accretion radius Racc^ such that particles that fall to 
smaller radii are assumed to be accreted by the BH and removed from the simulation. We discuss this in much more 
detail in Sec. 4.2. 

For a monatomic ideal gas, we can calculate the SPH internal energy as 



^1 



INT 



= Yl —Z^^i^iPi ^ = 9 X] rriiAipl 



2/3 



7-1 



Furthermore, the temperature is related to the specific internal energy Ui = ^Aip\ 

2 



2/3 



(10) 



via the ideal gas law £^int = ^ksT: 



3kB 



/^rrip 2/3 _ fJ^rrip Pi 
ks ks Pi 



(11) 



where rrip is the mass of a proton and p 
with mass fractions X = 0.7, F = 0.28. 



0.617 is the mean molecular weight, assuming that the disk is a plasma 



4. SPH SIMULATIONS 

In this section, we present dynamical calculations of kicked disks. The calculations of Sec. 4.1 treat the dynamics 
assuming the gas is optically thick, so that radiative cooling can be neglected, while the calculations of Sec. 4.2 assume 
that any heating that would have resulted from shocks is immediately radiated away. The dynamics of real disks will, 
depending on the system parameters, fall somewhere between these two extremes and likely harbor a rich dependence 
on position and time. 
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4.1. Calculations with Shock Heating 

To study the effect of the SMBH kick angle on the resulting disk evolution, we perform runs where the kick angles 
away from vertical are 15°, 30°, 45°, and 60°. We have also performed a number of test calculations to optimize 
various SPH-related parameters including the number of neighbors and the Courant factors (see Sec. 3), along with 
numerical convergence tests to guarantee the validity of our simulations and determine the parameters for our fiducial 
production runs. 

Our production runs are summarized in Table 1 below. The bound mass at the end of the simulation, M5 , is defined 
by Eq. 14. The angular momentum vector L5 of the bound matter is in the original coordinates, with the initial 
angular momentum of the disk in the z-direction. The L5 vector defines the ^''-direction used to construct cylindrical 
coordinates in the plots below. The tilt angle of the kicked disk is defined by the condition Ouu = diTccos{Lz;b/\Lb\)- 
The approximate maximum luminosity of the disk may be estimated from (^^f^) , though we note that we allow 
for the SPH particles to shock heat. 

Each run uses A^ = 5 x 10^ equal-mass SPH particles, and the number of neighbors is chosen to be 200 initially. All 
runs are started from the same relaxed disk configuration. To construct it, we lay down particles uniformly in space 
and use the local density as the basis for a Monte Carlo rejection technique. This configuration is relaxed for a time 
interval t = 160, during which we apply a drag force 

(Idrag = {v - VR)/trel (12) 

with trei = 0.8 as the chosen relaxation timescale and vr the exact rotation law satisfying Eq. 7 above, with vr = 

l{rc)/rc. 

TABLE 1 

Summary of the production runs performed. 



Kick angle 


n 


Bound mass 


l^bl 


Lb 


Tilt 


angle 


(°) 


Max. Luminosity 


15 




0.73 


0.543 


(0.135,0.009,0.526) 




14.3 




0.016 


30 




0.65 


0.444 


(0.166, 0.011, 0.412) 




21.9 




0.024 


45 




0.60 


0.359 


(0.146, 0.010, 0.327) 




24.2 




0.039 


60 




0.57 


0.296 


(0.107, 0.006, 0.276) 




21.3 




0.073 


None 




0.9999 


0.830 


a 




a 




0.015/0.004^^ 



^For the unkicked disk, |L^|, |L^| < 10"^, and thus Ouh < l.Oe - 6 as well. 

^For the unkicked run, there is a brief burst of internal energy generation when the dynamical effects are turned on, yielding an internal 
energy generation rate dEi^^/dt = 0.015, but thereafter the maximum steady state luminosity is dEi^^/dt = 0.004. 

Once the initial disk is relaxed, it is allowed to evolve dynamically until ikick = 0.8 before a kick is applied, except 
for a single unkicked control case we evolve to ensure that the physical effects we attribute to the kick are not merely 
an inevitable consequence of the dynamical evolution. In the discussion that follows, we define the time since the kick 
r^ as 

t^ ^ t tkick' \^'^) 

As shown in the bottom panel of Fig. 3, energy conservation is nearly exact for each of the runs, with total variation 
of no more than 0.03% in the total energy after the kick in any of the runs. Achieving this level of conservation is a 
consequence of two important components of the evolution scheme: the softened BH potential, described in Eq. A6, 
and the use of a Lagrangian-based variational scheme for evolving the smoothing length described in Sec. 3.2. The 
former, which may be justified given the finite spatial extent of an SPH particle, prevents particles on highly eccentric 
orbits that approach very closely to the BH from attaining spurious energy during the periapse passage. The latter, 
also used in Rossi et al. (2010), is required to allow for variable smoothing lengths without the energy varying on the 
same timescale as the smoothing lengths themselves. 

As can be seen in Fig. 3, the overall level of internal energy generation within the unkicked disk is approximately 
30% that of the most vertical kick configuration we consider (15° away from vertical), and roughly six times less that 
of the 60° kick simulation. We infer that while some of the heating we observe is an inevitable consequence of the disk 
evolution, the majority may be attributed to the kick and its aftermath, especially for cases where the kick is closer to 
the original disk plane. Similarly, the changes in the kinetic and potential energy seem to be almost entirely a result 
of the kick. 

By the end of our simulations, the kicked accretion disks have nearly reached a steady-state, as indicated by Fig. 3 
(for global quantities) and Fig. 4 (for bound matter). In general, the more oblique the kick, the more the resulting disk 
generates thermal energy, and, correspondingly, the deeper the potential energy well characterizing the disk. While the 
total kinetic energy is nearly uniform among all the kick angles we consider, we note that the bound mass is smaller for 
more in-plane kicks given the initial configuration we chose, and thus the specific kinetic energy of the disk increases 
with the obliquity of the kick. 

Fig. 4 shows the kinetic energy and relative mass of the bound portion of the disk after the kick for the different kick 
angles we considered as well as for the unkicked model. Clearly for the unkicked case, all the disk remains bound to 
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kick 



Fig. 3. — Evolution of the total potential energy (top panel), kinetic energy (second panel), internal energy (third panel), and total 
energy (bottom panel) for each of the runs. The kick occurs at f* = (see Eq. 13). Each panel shows the unkicked control run (violet 
dotted-dashed line), and kicked cases of 15° (red dotted line), 30° (green dashed line), 45° (blue long-dashed line), and 60° (solid black line) 
away from vertical. Energy preservation is almost perfect, with a total variation of no more than < 0.03% for any run, as shown in the 
inset plot in the bottom panel. 

the black hole, while for the different kick angles virtually all of the unbinding occurs at the moment of the kick. The 
exact bound fractions are determined by our choice of initial disk configuration; our bound disk masses are particularly 
sensitive to the angle subtended by the bound region of the disk at radii corresponding to the maximum surface density. 
Note that our definition of binding is the criterion 



^POT,i + ^KIN,i + ^INT,i 



rrii 



^. 



\viV ^Ap~ 



2/3' 



<o, 



(14) 



since the internal energy of an SPH particle on an otherwise bound trajectory would eventually lead to adiabatic 
expansion and unbinding of the constituent gas. Disk heating does lead to some additional unbinding of material in 
most runs, ranging from no additional mass loss at all up to 1.5% of the total mass, most of which occurs shortly after 
the kick. 

Shortly after the kick, at times < t* < 1, the net change in the disk's kinetic and potential energies is strongly 
dependent on the kick angle, as can be seen in Fig. 4. This is a purely geometric result that can be described in terms 
of a simple collisionless model (see Sec. 2.2): For a vertical kick, the entire mass of the disk is receding away from the 
BH immediately after the kick, and this remains nearly true for kicks near vertical. For more oblique kicks, where the 
bound component of the disk is drawn primarily from the part of the disk whose rotational velocity is aligned with the 
kick itself, a substantial part of the mass finds itself on orbits approaching the BH immediately after the kick before 
collisions and shocks begin to circularize the new orbits. Thus, for the oblique cases, we see an instantaneous decrease 
in the kinetic energy owing to the kick itself, followed by a rapid increase of larger magnitude as matter falls toward 
the BH. For more vertical cases, the effect is reversed, and we see an instantaneous jump in kinetic energy followed by 
a gradual turnaround and decrease. In both cases, the potential energy evolves accordingly, becoming more negative 
for the oblique cases relative to the vertical cases. 
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Fig. 4. — Evolution of the kinetic energy (top panel) and mass (bottom panel), relative to that of the initial disk, of the particles that 
remain bound to the SMBH. Conventions are as in Fig. 3, and the binding energy prescription is given by Eq. 14. 

It is well-understood from previous calculations that the post-merger disk will be substantially tilted with respect to 
the initial equatorial plane, so we chose a simple prescription to define the post-kick disk plane. Considering only the 
bound particles, as defined by Eq. 14, we calculate the angular momentum of the bound component of the disk with 
respect to the black hole, yielding the results shown in Table 1. Labeling this as the z-prime direction z\ we define x' 
to be the original x-direction with the parallel projection of z^ subtracted away, and the ^'-direction to be the cross 
product of the other two new coordinate directions: 



L, 



^/ -^bound 



X — {x - z'^z' 
\x — [x • z')z'\ ' 



y' = z' X X . 



In all of the plots that follow, radii are assumed to represent cylindrical radii in the primed coordinate system. 

To confirm the validity of the "gap-filling" model we discussed in Sec. 2.2, we show the SPH densities of the particles 
in the inner disk as a function of the cylindrical radius in Figs. 5 and 6. Turning first to the unkicked control model in 
the top panel of the former, we see that there is relatively little density evolution except at the innermost edge of the 
disk, where viscous dissipation of angular momentum leads to an accretion of particles toward the SMBH. A density 
peak does form at the center, but with smaller densities than the kicked runs at any given radius f < 0.1 during the 
duration of our simulations. 

Considering the kicked runs next, we observe that the filling of the gap proceeds as predicted above. A wave of 
particles at relatively high densities (p '^ 1 — 10) is observed moving inward at £* = 4 for the two most oblique kick 
angles we consider, namely 9 = 45° and 60°, particularly the latter. No such densities are ever found except in the 
immediate vicinity of the BH (f < 0.05) for the more vertical kicks {0 = 15° and 30°). Part of this effect is a simple 
matter of the larger post-kick periastron radii in the more vertical cases (see Fig. 1). Also, since the entire inner region 
of the disk remains bound in these cases, it forms an "inner ring" that will block any infalling matter from reaching 
radii f < 0.1 — 0.2. For the more oblique cases, not only does some matter accrete promptly, but a gap is formed in 
the inner region corresponding to the initially retrograde portion of the disk (i.e., the material initially located near 
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Fig. 5. — SPH particle-based density profile for the unkicked run (top panel), showing the maximuni (thick lines) and average SPH 
densities (thin lines) for particles binned with respect to cylindrical radius at £* = (red dotted lines), t* = 4 (green dashed lines), t* = 8 
(blue dot-dashed lines), and t* = 12 (black solid lines) . Note that the t* = configuration presents the common pre-kick state for all of 
our runs. There is a much more substantial density profile evolution in the kicked models (bottom panels). From left to right we show the 
evolution of our kicked models at t* = 4 (left), t* = 8 (center), and t* = 12 (right), with the maximum SPH density shown in the upper 
panels and the average SPH density in the lower panels. The curves follow the conventions of Fig. 3. 

(j) = —90°) that allows material to be channeled more easily toward radii f < 0.1. This resulting density enhancement 
surrounding the BH is quite persistent; at t^ = 12, which represents several orbital timescales for the innermost edge 
of the disk, there is a factor of five difference in the density at f < 0.1 between the 60° and 15° runs. We note that 
there is some "crystallization" and pattern formation of the particles located closest to the BH, which is inevitable 
when a small number of SPH particles are located near the edge of a density distribution, but that we still are able to 
resolve a smooth density field there. 
For the most oblique model, there is a clear oscillation in the maximum density at early times, which is highlighted 
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Fig. 6. — SPH particle-based densities for the innermost part of the accretion disks. All conventions follow those in Fig. 3 and the bottom 
panel of Fig. 5. 

in Fig. 6, showing the densities for the innermost region of the accretion disk. These features are an indication of 
the spiral density waves present in the disk and leave a clear imprint on the resulting thermodynamic evolution and 
emission properties we predict, as we discuss in more detail below. The general behavior of these density pulses, which 
move inward with time, is to increase the central density of the disk, as we can see by the closing of the gap between 
the densities for the 45° and 60° runs between t* = 8 and t^ = 12. For the more vertically kicked runs, where spiral 
waves are much weaker, there is very little sign of rapid accretion of material to the center. 

In setting up the runs shown here, it became evident that the BH softening potential we apply, Eq. A6, can play an 
important role in suppressing spurious energy fluctuations. Indeed, in simulations without a BH softening potential, 
we find that the innermost particles around the BH clump together, similarly to the so-called pairing instability of 
SPH (Price 2010) but with more than two particles per clump. These clumps form quasi-stationary "bubbles," where 
mutual pressure interactions keep any particle from approaching within about a smoothing length of the black hole. 
This behavior is robust against several different choices of the SPH smoothing kernel definition and evolution schemes 
for the smoothing length in time. When outside interactions finally "pop" the bubble, and other particles are able 
to flow inward to smaller radii, the measured kinetic and potential energies are seen to jump by substantial amounts 
because of a handful of particles, even though the total energy remained well conserved. 

For a slightly different view of the accretion disks, we show the azimuthally averaged radial surface density profiles 
I](f) in Fig. 7. We see that the surface densities are markedly different in the inner region, with the oblique kicks 
leading to persistently higher surface densities by at least a factor of five compared to more vertical kicks for r < 0.2 
throughout the course of the simulation. In the outer regions of the bound disk, the surface density trend is reversed 
but much less dramatic: the more vertical kicks have slightly larger surface densities at a given radius than the oblique 
kicks, but the variation is never more than a factor of two once the disk begins to relax again at t^ > 4. The ratio 
of the initial surface density to the postkick surface density is relatively constant over a wide range of radii, from 
0-4 ^ ^ ^ 1-5, but the postkick disk extends much further, since the same angular momentum exchange processes that 
funnel matter toward the BH at the inner edge of the disk also help to circularize it to larger radii at the outer edge 
of the bound region. 

Our results also allow us to make some rough conclusions about the opacities of our disks, though we note we do 
not include any radiation transport effects nor radiative cooling in the simulations of Sec. 4.1. Because the disks are 
hot and diffuse throughout, we expect the Thomson opacity for an ionized plasma to be a reasonable approximation. 
Whether or not the disk is optically thick (r >> 1) depends not only on the dimensionless surface density S but also 
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Fig. 7. — Radial profiles (averaged over azimuthal angle) of the surface densities (in log scale) for the initial disk model and the no-kick 
model (top panel) and for kick angles of 15, 30, 45 and 60 degrees; at t* = 4, 8, 12 (lower panels). 

on the values of the disk mass mdisk, the BH mass Mbh, and the kick velocity '^kick (see Eq. 3). For our reference model 
(^disk = 10^ M0, Mbh = lO^M©, and 'Ukick = lO^cm/s), the optical depth is slightly larger than unity throughout the 
pre-kick disk and slightly below throughout the post-kick disk, in which T^posti^) /^pre{^) ^ 0-2 — 0.4 for most of the 
area of the disk, 0.4 ^ r < 1.5. If the initial disk had a substantially smaller surface density, our model would predict 
that the post-kick disk would remain so as well, except in the very central region near the BH. Meanwhile, an optically 
thick initial disk should produce a slightly less thick disk after the kick, extending outward to nearly the edge of the 
bound component. 

While the azimuthal averages provide a clear picture of the global behavior of the disk, they do mask some of the 
more local phenomenon that develop after the kick. In Figs. 8 and 9 we show colormaps of the surface density at 
larger and smaller scales, respectively, with different color mappings between the two to allow for maximum contrast. 
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The large scale maps show the global expansion of the disk, though some care is required to interpret the results 
for the more oblique kicks. Most of the material that remains bound after the kick falls inward quickly to form the 
circularized inner disk extending out to f < 2.0. The lower arc that appears clearly for the 45° and 60° on the left 
(smaller x'-values) of the inner disk (and more subtly so for the 30° case) represents the unbound upper portion of 
the original disk (i.e., the material initially near (j) = 90°), with the density enhancements primarily a 2-dimensional 
projection effect. The material lies below (i.e., at lower ^''-values) that the bound component of the disk, except 
for the marginally unbound material where the stream connects to the inner disk. The upper arc visible clearly in 
the more oblique cases represents the marginally bound component of the disk from the upper edge of the original 
disk. This stream of material initially moves way from the BH after the kick while remaining level vertically with the 
inner disk before accreting toward the inner disk following a roughly ballistic trajectory. When it collides with the 
inner disk, it shock heats and circularizes, with an accretion rate that decreases gradually over time. In each of our 
runs, the shock fronts are never particularly sharp, certainly less so than the adiabatic 2-d thin-disk calculations in 
Corrales et al. (2010) which are themselves much more spread out than isothermal calculations in which there is no 
shock heating. Instead, because we allow the gas to heat as it shocks, the spiral patterns rapidly blur, leading to more 
extended density enhancements. The non-axisymmetry is strongest during the early phases of the simulation, and 
gradually fades as the disks relax and collisions circularize the fluid, so that by t* = 12 we see only minor deviations 
from axisymmetry, particularly near the outer edge of the bound region where infalling matter is still playing a role. 

In the smaller-scale surface density plots. Fig. 9, the role of the "gap" at f * = 4 is immediately apparent. For the 
60° kick, the center of the disk has already filled in, though the surface density is strongly non-axisymmetric even at 
very small radii. Significantly more matter is located at small separations for the 45° kick simulation, but a hollow is 
clearly visible around the BH. For the more vertical kicks, the gap is clearly present and very little matter is evident 
to begin filling it. By t* = 8, this influx of matter leads to a very sharp increase in the central density for the more 
oblique kicks that is not present in the more vertical ones. Finally, by t* = 12, the disk exhibits a much greater degree 
of axisymmetry, with only the most oblique kick case, in which the bound component of the disk is drawn almost 
entirely from one side of the pre-kick accretion disk, still retaining a marked angular dependence pattern. 

Turning to the thermal evolution of the disks, we see in the first panel of Fig. 10 the radially averaged thermal 
profile of our adiabatic initial model, for which T oc p^^^ . We find substantial heating very quickly in the center of the 
disk, with strong shocking throughout the region r < 1. The strong shock heating extends out through the majority 
of the bound component of the disk, except for the outermost regions at early times. Indeed, the sudden dropoffs in 
the temperature profiles corresponding to the more vertical kicks (15° and 30°) at f * = 4 occur toward the outer part 
of the bound region of the disk, not the unbound part of the disk. This is in accordance with our previous discussion, 
as the further reaches of the bound components for the more vertical kicks are initially sent outward in their orbits, 
and take longer to eventually collide with other fluid streams and shock. Over time, the temperature profile smooths 
out, and by £* = 12 we have essentially a single temperature profile, peaked in the center and then falling off more 
slowly at larger radii, that characterizes all of our kicked disk simulations. 

Based on the surface density and temperature profiles of our simulations (Figs. 7 and 10), it becomes clear that 
the differences in the global energies of the disks as a function of the kick angle are not the result of radically 
different temperature profiles, nor significantly different densities throughout the bulk of the disk, but rather from 
density differences in the innermost region of the disk. It is at small radii that all three energies take on their largest 
magnitudes, and the factor of 5 — 10 difference in the surface densities at these radii represent a significant fraction of the 
total energy, though not the total mass. Perhaps the clearest prediction from our simulations is that we expect oblique 
kicks to produce much more energetic signatures and substantially higher measured accretion rates for timescales of 
thousands of years for systems similar to our reference model and almost an order of magnitude longer for systems 
with kicks of roughly 500km/s. 

Our dynamical treatment assumes the gas is optically thick, so that the disk is allowed to shock heat without 
radiative cooling applied. Nevertheless, if we estimate the potential luminosity of the disk by assuming that the 
internal energy gains would be immediately radiated away but not affect the dynamics in any other way, we can 
convert the internal energy profiles from Fig. 3 into a luminosity by assuming that (^) , = ^^f^- The results are 
shown in Fig. 11, where we differenced over intervals 5t^ = 0.2 to minimize spurious noise. In all cases, we see an 
immediate but purely numerical luminosity peak when we end relaxation and begin the dynamical evolution. For the 
unkicked disk, we see only a very small, nearly constant luminosity over time. For the more vertical kicks, we find a 
small rise in luminosity followed by a gradual decline, but with very little temporal structure. For the more oblique 
kicks, we see both significantly larger luminosities as well as quasiperiodic emission spikes, especially for the 60° case, 
in which there are persistent oscillation amplitudes of tens of percent with a period of slightly longer than t = 1. Given 
the density patterns we observed, it seems clear that we are observing periodic shocking due to intersecting flows 
followed by accretion events as dense regions within the inner disk fall toward the BH while heating up significantly. 
The timescale roughly corresponds to the orbital timescale at the inner edge of the disk, with the strong m = 1 mode 
dependence of the fluid density (i.e., a single-arm spiral pattern) leading to increased shocking when the pattern wraps 
a full time around the BH. 

To confirm that these results are robust, and not merely a numerical artifact, we re-ran the 60° case using 1.25 x 10^, 
2.5 X 10^, and 10^ SPH particles in order to compare to our N = b x 10^ reference model. The luminosity we derive 
from each run is shown in Fig. 12. We see good agreement throughout the early phases, with slightly lower luminosities 
during the first peak for higher particle numbers because spurious numerical dissipation is slightly smaller. The drop 
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Fig. 8. — Projected surface densities (in log scale) for No-kick, and kick angles of 15, 30, 45 and 60 degrees (rows from top to bottom 
respectively); at £* = 4, 8, 12. 
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Fig. 9. — Projected surface densities (in log scale) for No-kick, and kick angles of 15, 30, 45 and 60 degrees (rows from top to bottom 
respectively); at £* = 4, 8, 12. 
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Fig. 10. — Radial profiles (averaged over azimuthal angle) of temperatures (in log scale) for the no-kick model (top panel), and kick 
angles of 15, 30, 45 and 60 degrees; at U = 0, 4, 8, 12 (from left to right). 



in luminosity at t* = 1 is a robust effect, as is the second luminosity peak at t* = 2. There is evidence for further 
oscillations in the luminosity, but not at a level we feel is conclusive. Still, given the strong convergence of the numerical 
results, it seems clear that the quasi-periodic luminosity is a physical prediction of our model, and not a numerical 
artifact owing to under-resolution of the accretion flows. 

4.2. Prompt radiation calculations and comparison to previous work 

Our assumption of Sec. 4.1 that the disk is allowed to shock heat differs from most previous works, particularly 
Corrales et al. (2010) and Rossi et al. (2010). We note that the assumption is closer to reality than it is often given 
credit for; indeed, while our default parameters from §2.1 indicates that the initial disk is optically thin, faster kicks, 
larger disk masses, or smaller black holes can potentially tip the balance toward optical thickness across much of the 
disk: see Eq. 3 and Fig. 7. In the innermost regions of the disk, where the majority of the luminosity actually arises, 
optically thick regions do develop, and it seems clear that immediate and passive cooling itself neglects important 
physical contributions to the dynamics of the disk. In order to bookend the likely observable emission from disks, we 
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Fig. 11. — Potential disk luminosities from our disk models, calculated as the time derivative of the internal energy as shown in Fig. 3. 
To smooth the data and minimize SPH discretization noise, we difference over intervals of 500 (No kick, 15° and 30°) or 1000 timesteps 
(45° and 60°). Besides the expected increase on average luminosity as the kick becomes more oblique, we also see substantial oscillations 
in luminosity for the more oblique kicks, especially the 60° case. 

have also run simulations using a "prompt radiation" formalism, in which the energy evolution equation, Eq. A2, is 
ignored, and we instead assume that all energy is immediately radiated according to Eq. A3, resulting in an adiabatic 
evolution. 

Calculations we perform using the prompt radiation assumption indicate that if the material does not get shock- 
heated, there are severe numerical challenges in evolving the inner region even when the BH potential is softened. 
Indeed, given the lack of thermal pressure support, we find generically that particles accumulate catastrophically 
around the BH, leading to a rapid slowdown of the evolution when the Courant timestep grows too small. Prior 
SPH works in which a prompt radiation assumption was used (see, e.g., Rossi et al. (2010)) have implemeted an 
accretion radius Racc^ inside of which all particles are immediately accreted and removed from the simulation. In 
Rossi et al. (2010), the accretion radius was set to the inner edge radius of the initial disk. Race = 0.1 (G. Lodato, 
private communication) . 

To allow a more direct comparison with Rossi et al. (2010), we have run two prompt radiation models, using a 60° 
kick , setting the accretion radius to be Race = 0.098, just inside the inner edge of our initial disk, and Race = 0.05, 
at half that radius. (We note that our angular convention is reversed from theirs: we refer to a vertical kick as ^ = 0° 
and an in-plane kick as 90°.) For such runs, strict energy conservation is impossible to maintain since the density 
contribution from accreted particles toward other particles is instantaneously removed in the middle of each timestep. 
Still, it is straightforward to calculate the implied luminosity rate, which we show in Fig. 13 for the two prompt 
radiation runs and our original run from Sec. 4.1 with shock heating included. 

Adding an accretion radius has several rather dramatic effects on the resulting luminosity. It smooths out the 
predicted luminosity, since high-energy interactions are suppressed, and increases the minimum timestep since particles 
never approach as closely to the BH. It also introduces an unphysical scale into the results, which becomes clearly 
evident when we compare results from the Race = 0.098 and Race = 0-05 cases. Indeed, we see an increase of over 50% in 
the predicted luminosity when the accretion radius is halved, since the characteristic energies of particles immediately 
before accretion, which represent a substantial component of the energy release, are roughly a factor of two greater 
based on dimensional arguments. As a result, conclusions about the potential luminosity of the disk unfortunately 
depend strongly upon the ad hoc choice of accretion radius. Also, we note this produces an unphysical accretion-driving 
process, since pressure support from the inner disk is removed resulting in strongly anisotropic pressure support for 
particles just outside the accretion radius. This leads to a persistent dependence in the long-term luminosity on the 
accretion radius, which is worrisome. 

A secondary effect of implementing an accretion radius is to eliminate periodicity in the luminosity of the disk for 
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Fig. 12. — Disk luminosities from our 60° cases with varying numbers of SPH particles, showing the results for N = 10^ (solid), 5 x 10^ 
(long-dashed; this is the same run shown in Fig. 11), 2.5 x 10^ (short-dashed) and 1.25 x 10^ (dotted) SPH particles. We see clear evidence 
for two peaks in the luminosity profle, at roughly t* = 0.7 and 2.1, and possible evidence for further modulation of the signal at later times. 



kicks of high inchnation angles, which we beheve results from excising from the computational domain the inner region 
of the disk where accretion flows with the highest energies typically interact. This is a problem for the 60 degree 
case in particular, as the post-kick periastron distances fall within the accretion radius for much of the initial bound 
component of the disk. 

We can now classify the effects of both shock-heating and inclusion of the innermost disk, as in our simulations 
of section 4.1, on the results of numerical simulations. Here, we see much more modulated emission produced from 
oblique kicks, as a result of high-eccentricity flows into the inner disk, while more vertical kicks yield much smoother 
luminosity profiles. In contrast, Rossi et al. (2010) found substantial oscillations for a vertical kick, which we suggest 
might result from modulations in the accretion rate into the inner disk but exclusion of the emission that would 
otherwise take place within f < 0.1. Our luminosities also peak at later times than their simulations indicated (see 
their Fig. 20). Among the possible explanations, the likeliest is that our simulations are most sensitive to colliding 
flows in the inner regions of the simulation, which require infalling material, whereas the simulations of Rossi et al. 
(2010) will have accreted much of the inner disk by the time such flows reach the center. As a secondary issue, the 
different initial disk configurations may also play a role, as we assumed an equilibrium model while theirs followed a 
power-law profile. The issue of adiabatic versus isothermal simulations also deserves mention, though spiral shocks are 
often stronger in isothermal models (Corrales et al. 2010). 

There is a precedent for the periodic luminosities we see, as a very similar profile appears in Fig. 4 of Shields & 
Bonning (2008), who studied a collisionless model for kicked disks that modeled emission via passive "collisions" 
assumed to occur when particles approached close to each other, though without any feedback on the resulting 
dynamics. They also found two strong luminosity peaks and a quasi-periodic luminosity profile afterwards, though 
for a 45° kick (we note that their disk model was different than ours, as it lacked pressure support, and that the 
normalization of time between their runs and ours differs by a factor of 2it). This strongly suggests that the kick 
luminosity in our simulations is driven by interacting flows, even in the presence of disk heating, and that these flows 
are not being modeled in simulations that excise the inner region of the disk. 

Although our initial disk model differs from that of Rossi et al. (2010), it still provides an opportunity to check 
our results alongside the "circularization" model proposed there, in which one calculates the approximate available 
energy for a fluid element within the disk by using its post-kick specific angular momentum / to infer a circularization 
radius Vc = L'^/GM and energy — GM/(2rc), and subtracting this away from the fluid element's initial energy after 
the kick (see their Eqs. 6 - 10). Rather than construct the energy estimate by a surface integral over the disk, we 
compute it particle-by-particle from the common disk configuration present immediately after relaxation. Our results 
are shown in Table 2, where the columns list the predicted "circularization energy" Ec^ the actual internal energy 
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Fig. 13. — Disk luminosities from our 60° case with the effects of shock heating included (solid curve), compared to two "prompt radiation" 
runs in which shock heating is ignored and luminosities calculated implicitly from Eq. A3. For the latter, we impose an accretion radius 
Race inside which all particles are assumed to be accreted by the BH and removed from the numerical simulation, choosing Race = 0.098 
(long-dashed, just inside the inner edge of the original pre-kick disk) and Race = 0.05 (short dashed, at half that separation). We see 
substantially larger luminosities for the prompt radiation cases, but with a strong dependence on the choice of the accretion radius, which 
is not a physically motivated parameter but rather a numerically chosen one. 

^INT at i^ = 12, and the difference between the internal energy of a kicked model and our unkicked reference model 
£'iNT;/cic/c = ^INT — ^iNT;no/cic/c cit t^ = 12, for cach of the kicked runs. We see that while the overall values are 
relatively close, indicating that the circularization model yields a good approximation for the order of magnitude of 
the energy that could be emitted by the disk, this very simple model does not yield quantitatively accurate predictions 
with respect to the kick angle dependence. Thus, while the formula is certainly useful for establishing the approximate 
disk luminosity given a reasonable timescale for emission, we do not see strong evidence that it can be extrapolated 
to angles that lie very close to the disk plane, for which the predicted energy rises like a power law in the out-of-plane 
angle (90° — 6k in our notation) to extremely large values (see Figs. 4 and 21 of Rossi et al. (2010)). 

TABLE 2 

Estimated energies available to the disk 



Kick angle (°) 


Circularization energy 


Internal energy at f* 


= 12 


Internal energy, corrected for kick 


15 


0.13 


0.12 




0.07 


30 


0.11 


0.16 




0.11 


45 


0.19 


0.24 




0.19 


60 


0.48 


0.32 




0.27 



5. DISCUSSION AND FUTURE WORK 

In this paper, we have studied the response of a quasi-equilibrium accretion disk around a SMBH that undergoes an 
impulsive kick, presumably because of a binary merger and the corresponding asymmetric emission of linear momentum 
in gravitational waves. While there are a number of sources that have been identified as candidate kicked disks with 
recoil velocities so large that they would be in the far tails of the kick velocity distribution, our results here are 
scale-free with regard to the kick, and may be applied to a broad swath of potential kicked systems. Indeed, while 
for our assumed reference model with Mbh = lO^M©, rridisk = lO^M©, and 'Ukick = 1000 km/s we find characteristic 
timescales of roughly 1000 years and disk luminosities and temperatures of up to lO^^erg/s and 10^~^K, respectively, 
the results should work for a wide range of masses and kick velocities. There is a trade-off, to be had, of course. 
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especially given the dependence of many of the quantities we investigate on the kick velocity. Indeed, cutting the kick 
velocity by half and leaving the masses of the BH and disk unchanged increases the timescales we consider by a factor 
of eight while cutting the energies and temperatures by a factor of four and thus the luminosities by 32, indicating 
that there should be a definite observational bias toward the larger kicks. 

The code we introduce uses 3-dimensional SPH techniques, and is modified to incorporate a Lagrangian ("grad-h") 
prescription that allows energy to be conserved to high precision over the course of all our runs. We also introduce 
a Lagrangian-based black hole smoothing potential, which proves critical in allowing us to avoid numerical issues 
associated with point-like potentials. One of the most important purely numerical conclusions of this work is that 
black holes must be handled extremely carefully in SPH, since they can introduce particle clustering instabilities that 
are both conservative but highly unphysical. This treatment, along with the inclusion of shock-heating, has allowed 
us to perform what we believe to be the first 3-d SPH simulations of these systems that do not require imposing an 
ad hoc "accretion radius" inside of which SPH particles are removed from the simulation. 

In order to examine the phase space of post-kick disks, we have varied the kick angle of the SMBH with respect to 
the initial disk orientation, which has a large effect on the resulting evolution. More oblique kicks, i.e., those most 
oriented toward the equatorial plane of the original disk, produce substantially higher peak luminosities for a given BH 
and disk mass and kick velocity, with roughly a factor of four gain between kicks oriented 15° away from vertical and 
those oriented at 60° from vertical. Assuming an astrophysical context in which the disk is aligned with the SMBH 
binary orbit prior to merger, it is unlikely that more oblique kicks will actually yield more luminous events, however. 
Kick velocities are systematically higher for more vertical kicks (Lousto et al. 2010), and given the vl-^j^ dependence 
of the luminosity, the brightest kicked disks are likely to be those with the highest speed kicks, with the kick angle 
playing a secondary role. 

We find more rapid luminosity peaks appearing for more oblique kicks, which is the opposite result from a previous 
set of simulations that considered power-law density profile disks (Rossi et al. 2010). Based on our 60° calculation, we 
attribute this to the rapid filling of the innermost region of the disk by material kicked into high eccentricity orbits 
with extremely small periastron distances (see Fig. 1), which fiows inward from the inner edge of the pre-merger disk; 
we note that this pattern differs from that found in Fig. 15 of Rossi et al. (2010), likely due to the combined effects 
of shock heating and their excision of particles from the innermost regions of the disk. Clearly, one of the important 
conclusions of this work and others is that fiows of material toward the SMBH after the kick can release tremendous 
amounts of energy, but need to be modeled very carefully to derive reasonable light curves and spectra, pushing the 
limits of current numerical simulations. 

The surface density oscillations and the resulting quasi-periodic emission we observe in time for the most oblique 
kicks will make an interesting topic for further research. Based on our "gap-filling" model it seems clear the dynamics 
of the inner disk play an important role in the dynamical evolution and EM emission from the disk, and it would be 
highly useful to extend the same simulations on both sides of the kick to explore the full history of a circumbinary/post- 
kick disk. In particular, it will be important to accurately resolve the inner edge of the circumbinary disk, where the 
tidal field of the binary should have swept out a gap (Schnittman & Krolik 2008), and to understand how infalling 
material from that region interacts with the more tenuous tidal streams of matter in the process of accreting onto the 
SMBH prior to their merger (Macfadyen & Milosavljevic 2008; Hayasaki et al. 2007; Hayasaki et al. 2008). It will also 
be useful to incorporate a more thorough treatment of the radiative evolution of the disk, since a proper treatment 
of radiative cooling and disk opacities will break the scale-invariance of the results and allow for much more accurate 
predictions of observable phenomena. 
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velocity relationship, and Fabio Antonini, Manuela Campanelli, Evgenii Gaburov, Julian Krolik, Carlos Lousto, David 
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NCSA and TACC under grant numbers TG-PHY060027N and TG-AST100048 and has made use of both the SPLASH 
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APPENDIX 

SPH EVOLUTION SCHEME 

Our SPH evolution scheme uses a variable smoothing length approach following the formalism described in Springel 
& Hernquist (2002) and Monaghan & Price (2001). Defining a particle-based parameter 
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while our evolution equation for the entropic variable Ai is 

For the "prompt radiation" simulations discussed in Sec. 4.2, we hold Ai fixed and ignore Eq. A2, though we do evolve 
Eq. AI as written. Energy losses are calculated by integrating the implied change in the entropic variable over time 
^^, and inserting this into the expression 

^^^^ = E -^MdA/dtU,HedPr' = E E ^^^:>^i^i - V.) • "^iWiAhi) . (A3) 

i i j 

In the above, Wij{hi) is the normalized second-order accurate kernel function introduced in Monaghan & Lattanzio 
(1985), used widely throughout the SPH community, and /i^, p^, rrii, and Pi are the SPH particle smoothing lengths, 
densities, masses, and pressures, respectively. The artificial viscosity term Uij is discussed below. We assume a 
smoothing length-density relation in the form 

h^=(-r^+hpy') (A4) 

\ ''"max / 

where hmax = 50. The bi values are chosen so that each particle should have ~ 200 neighbors given the initial density 
profile of the disk and are updated to maintain this condition during relaxation. Once the dynamical phase of the 
evolution begins, we hold bi fixed and solve implicitly at each time step for the proper smoothing length and density 
that satisfy equation (A4). 

Because the self-gravity of the disk is ignored here, the only gravitational force acting on the particles comes from 
the black hole. We assume instantaneous Newtonian gravitation, neglecting retardation effects from the moving black 
hole since our characteristic speeds are small fractions of the speed of light. Although we treat the black hole as a 
pure point mass without any intrinsic softening, the gravitational force on any SPH particle within two smoothing 
lengths of the black hole is softened according to the mass distribution of that SPH particle itself, as described by its 
smoothing kernel. To do so, we follow the formalism of Price & Monaghan (2007), who use a variational approach to 
derive equations of motion that properly account for variable smoothing lengths. In particular, we start by writing 
the gravitational part of the Lagrangian as 

^grav = - ^ rui^i = -GMbyl ^ miifi{hi). (A5) 

i i 

The last equal sign in equation (A5) defines the gravitational potential ^i of particle i. Here (pi{hi) refers to (p{ri^ hi) 
where r^ = Ir^l is the distance of particle i from the BH and 

0<g<l; 

^{r.h) = { (-l + ik + l^'-^' + ^^'-^^O/^' l<q<2; , (A6) 

q>2 

with q = r//i, is the gravitational potential associated with the usual SPH smoothing kernel (e.g. Price & Monaghan 
2007). 
Following the approach of §3 of Price & Monaghan (2007), but with our Lagrangian, we find 
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-h^' 


+ hi") /h 
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+ h'- 


- l' + i-ol' 



^(grav) ^ _GMBHVi(pi(/l,) - GMbh Y1 ' 






(A7) 



The first term on the right hand side of equation (A7) is the usual softened gravitational acceleration, while the 
remaining terms allow for variable smoothing lengths and preserve energy conservation. We note that one of these 
correction terms vanishes when the kernel of the relevant SPH particle does not overlap with the BH, because, given 
equation (A6), d(p/dh = whenever r > 2h. Summing equations (AI) and (A7), the total acceleration of particle i is 
calculated as 

Vi = -GMBn^^ip^(h,) - ^ mj [T^ V, W,,- (hi) + T,-,V, W^ (hj)] , (A8) 

j 



where 



^ij GMbh dcpi dhi 

"'^ ~ l^ipf ^ 2 ^i dhi d Pi 



and where, given equation (A4), we use dhi /dpi = —bih'^/{3p/ ). 
We implement the AV form 




^iJ = ^^+^] i-o'Pij + ppI) fi , (Aio) 



Disks around kicked black holes 



23 



with a = P = 1. Here, 



f^ij 



{Vi-Vj)-{Ti-Tj) 



hij{\ri-rj 
0, 



IV/^^ 



T]^ 



, if {vi -Vj) • (r^ -r^-) < 0; 

if{Vi-Vj)-{Ti-Tj) >0, 



(All) 



with 77^ = 10 ^, and the Balsara switch fi for particle i is defined by 



/. 



|V-v|, + |Vxv|, +r?'ci//i, ' 



(A12) 



^ preventing numerical divergences (Balsara 1995). The function fi approaches unity in regions of strong 
|V-v|i>>|Vxv|i) and vanishes in regions of large vorticity (|Vxv|i>>|V-v|i). Consequently, our 



with V = 10 
compression 

evolution equations have the advantage that the artificial viscosity (AV) is suppressed in shear layers. We note that 
the AV term is not symmetric under interchange of the indices i and j (that is, 11^^ ^ n^^), because the switch fi is not 
symmetrized in equation (AlO). Such an approach reduces the number of arrays shared among parallel processes. As 
the term in square brackets in equation (A8) is antisymmetric under the interchange of particles i and j, momentum 
is clearly conserved in every interaction pair. Similarly, it is straightforward to show total energy is conserved by our 
evolution equations: ^- rriiiyi • v^ + d^i/dt + dui/dt) = 0, where the specific internal energy Ui = AipJ~ /{j — 1). 

The evolution equations are integrated using a second-order explicit leap-frog scheme. For stability, the timestep 
must satisfy a Courant-like condition. Specifically, we calculate the timestep as 



For each SPH particle i, we use 



and 



At = Mmi{Ati^i,At2,i). 
hi 



AUi = C 



NA 



MdiXj (TijPif^'^ 



At2,i = Cn,2 



hi 



1/2 



(A13) 
(A14) 

(A15) 



The Maxj function in 



For the simulations presented in this paper, the Courant factors Cn,i = 0.4 and Cn,2 = 0.05. 

equation (A 14) refers to the maximum of the value of its expression for all SPH particles j that are neighbors with i 

The denominator of equation (A14) is an approximate upper limit to the signal propagation speed near particle i. 



SPH EXPRESSION FOR THE SURFACE DENSITY 



We adopt a kernel function 



W{r, h) = 



tt/i^ 



0, 



a)' + !a) 



where we have defined 
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The function Wout is just the kernel function in the outer regions of the compact support (1 < r/h < 2), while Win is 
the difference between the kernel function in the inner and outer domains. The distance r from the center o f an SPH 
particle to the points on a line passing vertically with horizontal offset p from that center is given by r = ^/Z^ 

where Z is the vertical offset. Thus, we may define Zout = \/4/i^ — p^, and if p < h the quantity Zin = \/h'^ — p 
define the integration bounds for the kernel-based surface density 



2 to 



^=E 



rrii 



W(r, hi)dZ. 



(Bl) 



In particular, a particle i that is passed through in only the outer part of its kernel by the line being integrated along 
(so hi < p < 2hi) will contribute 



^out,i — 



I/.? 



^'» (if) -="'(¥) 



'J . f ^out 

2^H^ 



1 . f ^out 



(B2) 
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to the surface density. While if p < hi then particle i contributes 



^in,i — ■^o 



TTh^ 



Here we make use of the vertical symmetry of the kernel function to define the following integrals: 
h{z) = hio (^) = / dZ = 2z 

hiz) = hHi (^^^ = f r{Z)dZ=j ^JZ^ +p^ dZ=^ (Z\/Z-^ + p^ + p'' In \z + ^ Z^ + /j2J ) 



(B3) 



z\/ z'^ ^ p^ ^ (? In 



Z + yj Z^ + p^ 



{Z'^p')dZ=^-^{z'^3p') 



h{z) = hH,Q=j r'dZ- 



In 



:+y^^T7 
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